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It is known that solutions of Richardson equations can be represented as stationary points of 
the "energy" of classical free charges on the plane. We suggest to consider "probabilities" of the 
system of charges to occupy certain states in the configurational space at the effective temperature 
given by the interaction constant, which goes to zero in the thermodynamical limit. It is quite 
remarkable that the expression of "probability" has similarities with the square of Laughlin wave 
function. Next, we introduce the "partition function", from which the ground state energy of 
the initial quantum-mechanical system can be determined. The "partition function" is given by a 
multidimensional integral, which is similar to Selberg integrals appearing in conformal field theory 
and random-matrix models. As a first application of this approach, we consider a system with the 
constant density of energy states at arbitrary filling of the energy interval, where potential acts. In 
this case, the "partition function" is rather easily evaluated using properties of the Vandermonde 
matrix. Our approach thus yields quite simple and short way to find the ground state energy, which 
is shown to be described by a single expression all over from the dilute to the dense regime of pairs. 
It also provides additional insights into the physics of Cooper-paired states. 

PACS numbers: 74.20.Fg, 03.75.Hh, 67.85. Jk 

I. INTRODUCTION 

Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity- is based on the reduced interaction potential, which 
couples only electrons with opposite spins and zero total momentum, while the interaction amplitude is taken to be 
momentum-independent. It is known for a long time that Hamiltonians of this kind are exactly solvable^^— . Namely, 
they lead to Richardson equations, which are nonlinear algebraic equations for the set of complex numbers Rj{j = l,-.-, 
N), where TV is the number of Cooper pairs in the system; the energy of the system is given by the sum of all R's. 
Originally, Richardson equations have been derived by solving directly the Schrodinger equation. However, they can 
be also recovered through an algebraic Bethe-ansatz approach^. The resolution of the system of Richardson equations, 
in general case, is a hard task^. These equations are now intensively used to study numerically superconducting state 
in nanosized superconducting systema^^— . They are also applied in the nuclear physics, see e.g. Ref.- *^ 

Recently, Richardson equations have been used^i to explore connections between two famous problems: the one-pair 
problem solved by Cooper— and the many-pair BCS theory of superconductivitj*i. The essential ingredients of the 
Cooper model are the Fermi sea of nonintcracting electrons and the layer above their Fermi energy, where an attraction 
between electrons with up and down spins acts. Two additional electrons are placed into this layer. It is also assumed 
that the energy density of one-electronic states within this potential layer, which has a width equal to the Debye 
frequency, is constant. Cooper then was able to solve the Schrodinger equation for two electrons. He found that the 
attraction, no matter how weak, leads to the appearance of the bound state. Cooper model was an important step 
towards a microscopic understanding of superconductivity. In contrast, BCS theory of superconductivity considers a 
many-pair configuration, which includes the potential layer half-filled. Traditional methods to tackle this problem is 
either to use a variational approach for the wave function^ or to apply Bogoliubov canonical transformations^*. In Ref. 
fii] , an arbitrary filling of the potential layer has been considered. Obviously, by increasing a filling, one can attain a 
many-particle configuration starting from the one-pair problem. Although such a procedure seems to be unrealistic 
from the view point of current experimental facilities, it allows one to deeper understand the role of the Pauli exclusion 
principle in Cooper-paired states, as well as to analyze a correspondence between the single-pair problem and BCS 



2 



condensate. This gedanken experiment can be also considered as a toy model for the density-induced crossover 
between individual fermionic molecules, which correspond to Cooper problem, and strongly overlapping pairs in BCS 
configuration^^"^^ . 

In order to find the ground state energy for TV pairs, in Ref. fii], a new method for the analytical solution of 
Richardson equations in the thermodynamical limit was proposed. Rigorously speaking, this method is applicable to 
the dilute regime of pairs only, since it assumes that all Richardson energy-like quantities are close to the single-pair 
energy found by Cooper. By keeping only the lowest relevant correction to this energy, the expression of the ground 
state energy was found to read as 



where epo is the Fermi energy of noninteracting electrons (or a lower cutoff), f2 is the potential layer width (which 
corresponds to Debye frequency in Cooper model), a — exp{—2/pV), p being a constant density of energy states and 
V an interaction amplitude. An analysis of this equation shows that first two terms in its right-hand side (RHS) 
correspond just to the bare kinetic energy increase due to N pairs, while the third term gives the condensation energy. 
This condensation energy has a quite remarkable structure: it is proportional both to the number of pairs N and to 
the number of empty states pVL -|- 1 — iV in the potential layer. Moreover, in the limit N ^ \^ the condensation energy 
per pair is exactly equal to the binding energy of a single pair, found by Cooperl^. Thus, a pair "binding energy" in a 
many-pair configuration appears to be simply a fraction of the single-pair binding energy, the reduction being linear 
in N and proportional to the number of occupied states. 

Although Eq. (1) was derived for the dilute regime of pairs, it is in a full agreement with the mean-field BCS result 
for the ground state energy (originally derived for the half- filling configuration). BCS approach has been also applied 
to the case of arbitrary filling of the potential layer— li^, this analysis being similar to BCS-like models developed in 
Refsji^ii^. The result for an arbitrary filling was also shownii^i^ to be consistent with Eq. (1). First two terms in 
the RHS of Eq. (1) must be identified with the normal state energy, as it appears in BCS theory. The third term is 
exactly the condensation energy, apart from the fact that in BCS theory it is traditionally written in terms of the gap 
A that masks a link with the single-pair binding energy. A detailed discussion of these issues can be found in Refs. 



In the recent paper [IS], the dilute-regime approach to the solution of Richardson equations^^ has been advanced 
to take into account the next-order term of the expansion, which is needed in order to get rid of the dilute regime. 
It was found that the corresponding contribution to the ground state energy is underextensive, i.e., negligible in 
the thermodynamical limit. This means that Eq. (1) derived in the dilute regime of pairs is most likely to be 
universal, which implies that mean-field BCS results are exact in the thermodynamical limit, in agreement with 
earlier conclusions''^'^^^. Finally, very recently the dilute-limit procedure has been extended by taking into account 
all the terms in the expansion^^ . It was concluded^^ that corrections beyond the dilute-limit result of Ref. are 
indeed underextensive. 

The aim of the present paper is to find a ground state energy using Richardson equations without utilizing any 
asymptotic expansions from the single-pair configuration and to address a validity of Eq. (1) (beyond both a mean-field 
approximation and dilute-limit expansions of Richardson equations). For that, we suggest another way to evaluate 
Richardson equations analytically in the thermodynamical limit. We start with the so-called exact electrostatic 
mapping: It was noted p] that Richardson equations can be obtained from the condition of a stationarity for the 
"energy" of the two-dimensional system of charges, logarithmically interacting with each other, as well as with 
the homogeneous external electric field. In the present work, we suggest to make one more step and to consider 
"probabilities" to find a system of charges in a certain configuration at the effective temperature, which is just equal 
to the interaction constant V . This constant goes to zero in the large-sample limit, such that the "energy" landscape is 
very sharp in the vicinity of its stationary points. The logarithmic character of the particle-particle interaction energy 
in two dimensions together with our choice for the effective "temperature" leads to the rather compact expression 
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of "probability". Note that it has similarities with the square of Laughlin wave function appearing in the theory of 
fractional quantum Hall effeclj^^. Next, we introduce a "partition function" and identify the ground state energy of 
the initial quantum-mechanical system with the logarithmic derivative of this function with respect to the inverse 
effective temperature. This "partition function" is actually represented by a multidimensional coupled integral. 
Coupling between integration variables turns out to be similar with that for Sclberg integrals (Coulomb integrals) 
appearing in conformal field theory (see, e.g., Ref. [^^]), in certain random matrix models (Dyson gas^S.)^ and in 
growth problems^^<. At the same time, the structure of the integrand for each variable indicates that it can be also 
classified as the Norlund-Rice integral. 

A connection between Richardson equations and conformal field theory has already been addressed in Ref. [^^], 
where it was shown that BCS model can be considered as a limiting case of the Chern-Simons theory. Indeed, it is 
easily seen that the structure of Richardson equations has similarities with the structure of Knizhnik-Zamolodchikov 
equationsSSi. At the same time, Chern-Simons theory also plays a very important role in the description of the 
quantum Hall effect ^] . That is why there are similarities between the expression of " probability" , as introduced in 
the present paper, and Laughlin wave function. Note that conformal field theories attract a lot of attention, since 
it is assumed that they can provide a unification between quantum mechanics and the theory of gravitation (anti de 
Sitter /conformal field theory correspondence). 

In the present paper, we restrict ourselves to the case of a constant energy density of one-electronic states. We 
found an efficient way to calculate the "partition function" by transforming the integral to the sum of binomial type. 
Such a transformation is possible due to the Norlund-Rice structure of the integrand. The sum is evaluated using 
properties of the determinant of the Vandermonde matrix, this determinant being responsible for the coupling between 
summation variables. We then are able to obtain Eq. (1) and to prove its validity. We also suggest a "probabilistic" 
qualitative understanding of this result, which is related to the factorizable form of "probability". We found that the 
"probability" for the system of pairs, feeling each other through the Pauli exclusion principle, can be represented as 
a linear combination of products of "probabilities" for single pairs, each pair being in its own environment with the 
part of the one-electronic levels absent. These missing levels form two bands in the bottom as well as at the top of the 
potential layer for each pair, such that the sum of energies of single pairs for each term of a factorized "probability" 
is the same. 

Another method to evaluate the multidimensional integral, used in the present paper, is to integrate it through the 
saddle point corresponding to the single pair. This approach is actually based on a well-known trick used to compute 
binomial sums by transforming them to Norlund-Rice integrals, which can be tackled by a saddle-point method^. In 
our case, this procedure appears to be quite similar to the solution of Richardson equations in the dilute regime of 
pairs, as done in Ref. fii]. It gives rise to the expansion of energy density in powers of pairs density. However, we were 
unable to perform calculations along this line beyond few initial terms due to the increasing technical complexity of 
the procedure. Nevertheless, we show that only the first correction to the energy of noninteracting pairs is extensive, 
while two others are underextensive, this condition being necessary for Eq. (1) to stay valid. Thus, the first method, 
based on manipulations with binomial sums, turns out to be much more powerful and, moreover, technically simpler. 
We believe that the method, suggested in the present paper, is applicable to situations with nonconstant density of 
energy states, as well as to other integrable pairing Hamiltonians, for which electrostatic analogy exists^i^^. 

We also note that there already exists a method for the analytical solution of Richardson equations in the ther- 
modynamical limit^i^i^. For the case of constant density of states, this method, up to now, has been applied to 
half-filled configuration only, for which its results agree with BCS theory. In Refj^, it was also used for any pair 
density and for the dispersion of a three-dimensional system. We here consider arbitrary fillings at constant density 
of states. In addition, the method of Ref. f^] assumes that energy-like quantities, in the ground state, are organized 
into a one-dimensional structure in the complex plane. This assumption is based on the results of numerical solutions 
of Richardson equations^. In this paper, which is purely analytical, we would like to avoid using this assumption, 
that is why we have constructed another method, which, moreover, is technically simple and provides some additional 
insights to the physics of Cooper pairs. Nevertheless, results of both approaches for the ground state energy do 
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coincide with each other, as weh as with results of the mean-field BCS theory. 
This paper is organized as follows. 

In Section II, we formulate our problem and we introduce a basis of our "probabilistic" approach. 

In Section III, we find the ground state energy in a rather simple way, by using a representation of the "partition 
function" through the coupled binomial sum. We also suggest "probabilistic" interpretation of the obtained result. 

In Section IV, we tackle integral entering the "partition function" by integration through a single-pair saddle point 
and we also establish a connection with the approach of Refs. fiiii^]. 

We conclude in Section V. 

II. GENERAL FORMULATION 
A. Hamiltonian 

We consider a system of fermions with up and down spins. They attract each other through the usual BCS reduced 
potential, coupling only fermions with zero total momenta as 

V = -VY^ ak/^alk'^a-kiakt- (2) 

k,k' 

The total Hamiltonian reads a.s H = Hq + V , where 

^^0 = ^ ^k (aj^^flkt + a-l^<^kij ■ (3) 

k 

It is postulated that the potential V acts only for the states with kinetic energies Ek and Ek' located in the energy 
shell between efq and + il. In BCS theory, the lower cutoff epo corresponds to the Fermi energy of noninteracting 
electrons, while is Debye frequency. We also assume a constant density of energy states p inside this layer, which is 
a characteristic feature of a two-dimensional system. For a three-dimensional system, this is justified, provided that 
n Epg. Thus, the total number of states with up or down spins in the potential layer is Nq = pH, while these states 
are located equidistantly, such that £k = £Fo+ Ck, where ^k runs over 0, 1/p, 2/p, fl. Energy layer accommodates 
N < Nn electrons with up spins and the same number of electrons with down spins. These electrons do interact via 
the potential given by Eq. ([2]). 

In this paper, we restrict ourselves to the thermodynamical limit, i.e., to A — oo, where A is the system volume. 
In this case, p ~ A and V ~ A~^, so that the dimensionless interaction constant, defined as v — pV, is volume- 
independent, ~ A°. The same is valid for epg and fl: hence Nn ^ A. The number of Cooper pairs scales as TV ~ A; 
consequently, filling N/Nq is volume- independent, ^ A°. We treat arbitrary fillings of the potential layer N/Nq, 
while the traditional BCS theory deals with the half-filling configuration. Studies of arbitrary fillings help one to 
reveal an important underlying physics, which is not easy to see when concentrating on the half-filling configuration, 
which is quite specific. 



B. Richardson equations 

It was shown by Richardson that the Hamiltonian, defined in Eqs. ([2]) and ([3]), is exactly solvable. The energy of 
N pairs is given by the sum of N energy-like complex quantities Rj (j = 1,-.., N) 

N 

En = Y.R,. (4) 
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These quantities satisfy the system of N coupled nonhnear algebraic equations, called Richardson equations. The 
equation for Rj reads as 

where the summation in the first term of the RHS of the above equation is performed for located in the energy 
interval, where the potential acts. Note that the number of pairs enters to the formalism through the number of 
equations that is rather unusual. The case TV = 1 at A oo, corresponds to the one-pair problem solved by Cooper. 
The fully analytical resolution of Richardson equations, in general case, stays an open problem. 

C. Electrostatic analogy 

Let us consider the function Eciassi{Rj}), given by 

i^cia..({i?j}) =2 K^i?ei?,+F^ln|2ek-i?j|-2t/ ^ \n\Ri - R,\\ . (6) 
This function can be rewritten as 

E.iass{{Ri]) - W{{R,-\) + W^({i?*}), (7) 

where 

= E + ^ E (2ek - R,) - 2y ^ In {Ri - R,) . (8) 

j j,k 3,l,i<l 

Richardson equations can be formally written^ as stationary conditions for W{{Rj}): dW{{Rj})/dRj = 0. It is easy 
to see that Eciass{{Rj}) represents an energy of N free classical particles with electrical charges located on the 
plane with coordinates given by (Re Im Rj). These particles are subjected into an external uniform force directed 
along the axis of abscissa with the strength —2. They are attracted to A^a fixed particles each having a charge —\/V 
and located at ek's. Free charges repeal each other. Richardson equations are equivalent to the equilibrium condition 
for the system of N free charges. 

D. "Probabilistic" approach 

A key idea of the approach we here suggest is to switch from the "energy" of the system of charges to the "proba- 
bility" S{\RjY) for this system to be in a certain state at effective "temperature" Tejf given by the simple condition 

fcsTe// ^ 1/, (9) 

where 

5(R}) = .xp(-l!2ffll)). ,10, 

Taking into account Eq. for W(\Rj\), we obtain for a nicely compact expression 



nf=ink(2ek-i?,)'''n ^ 



sm^) - ^r^": „ exp , (11) 



which has obvious similarities with the square of Laughlin wave function at filling 1, due to the factor J^^ ^ ^^^{Ri—Rj) 
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In principle, it could seem more reasonable to use Eciass{{Rj}) instead of W^({i?j}) in the definition of "probability" . 
Indeed, the corresponding function exp {—Eciass{{Rj}) /V) is real-valued and positive, so that by its properties it is 
closer to a usual distribution function, as compared to S{{Rj}). However, such a function is not meromorphic, so 
that Cauchy theorem is not applicable; therefore, it is not so useful for the reasons, which will be clarified below. 

The important fact is that the effective "temperature" T^ff goes to zero in the thermodynamical limit as ksT^ff ^ 
A~^. This makes W{{Rj})/V extremely large by its absolute value, while the landscape oi S{{Rj}) is very sharp in the 
vicinity of stationary points of W{{Rj}). Hence, it looks attractive to try extracting an information about stationary 
points of W{{Rj}) by using integration techniques. We briefly illustrate this idea. Let g{x) be a function of the variable 
X, which has a sharp maximum at x ^ Xq. It is possible to flnd xq approximately without solving equation g'{xo) = 
explicitly, but through the integration. Namely, we consider a ratio of integrals J x exp((7(x))(ia;/ J exp{g{x))dx with 
Xq located "deeply enough" inside the integration interval. The dominant contribution to both integrals is provided 
by a neighborhood of xq, so that we expect that their ratio is close to Xq- In particular, if a second and higher-order 
derivatives of g{x) at xq are proportional to some large parameter, for example to A, then it is easy to show, by 
performing Taylor expansion of g{x) around xq, that the error in the determination of xq through the above ratio is 
of the order of A^^. 

However, it is easily seen that stationary points do not necessarily correspond to minima of "energy", but they 
rather give its saddle points, so that equilibrium positions of free charges are not stable. The easiest way of seeing 
it is to consider a one-pair problem, for which equilibrium R is located on the real axis, as indicated in Fig. 1. The 
energy of the system Eciass{R) (as well as the real part of W{R)) increases, if we start to move R out of equilibrium 
from the real axis in a perpendicular direction, but it decreases upon motion along the real axis in both directions 
(see Section IV for more details). Hence, to apply a saddle-point method, we should use an integration path passing 
through the stationary point, as shown in Fig. 1 by line 1. This procedure could be seen as useless from the viewpoint 
of determination of a saddle point position, since we already need to know it to apply this technique. However, this 
is not true, because we use W{{Rj}) instead of Eciass{{Rj}) in the definition of the "probability" S. Function S is 
meromorphic, which means that the result of integration is independent on an integration path for all paths which 
can be continuously transformed to each other without crossing any pole of S. Hence, we can use a variety of paths, 
which start at —ioo and end at +ioo, where S vanishes. By using these "nonlocal" properties of S, we can reconstruct 
an information on the unknown location of a saddle point in the complex plane. Therefore, we are led to consider the 
ratio of integrals defined as 

jRS{R)dR 
^= jSiR)dR ' 

where integration is performed in a complex plane from —ioo to +ioo. From one point of view, by passing integration 
path through the saddle point, we must recover with accuracy 1/p the position of the saddle point. From another 
point of view, we may use any integration path, provided it avoids all poles in the same way as line 1 of Fig. 1. Taking 
into account Eq. (Ilip for S{{Rj}), we may rewrite Eq. (|12p in an equivalent form as 

i^^-^mz, (13) 



where 



For a single-pair problem, S reads as 



Z = J S{R)dR. (14) 



exp{-R/V) 



Figure 1: The schematic plot of complex values of R relevant for a single-pair problem. Filled circles show locations of one- 
electronic levels. Open circle indicates an equilibrium position of R for the ground state. Line 1 shows an integration path 
through this point in the direction of the steepest descent of —Re W . Lines 2 and 3 represent examples of integration paths, 
for which the results of integration will be the same as for the line 1. Dashed line depicts the auxiliary semicircle introduced 
to use a residue theorem. 



We then perform a partial- fraction decomposition of 1/ nn=o(^^-Po — R + 2n/p) by using standard rules and rewrite 
S as 



where C is an irrelevant TZ-independent constant (C — (p/2) " (iVfj!) which will be dropped below, as well as 
similar irrelevant factors, while 



is a binomial coefficient. We then substitute Eq. (IT6l) to Eq. p4|) and perform an integration. This can be done by 
various methods. For instance, we may perform integration in a straightforward way, for the integration path along a 
straight line, like line 2 in Fig. 1. We then see that, of course, the result of the integration depends only on how many 
poles of function S, as given by Eq. (jl6p (corresponding to the positions of one-electronic levels 2£k in the complex 
plane), is at the right (or left) side of the integration path and it stays the same as long as the number of poles at the 
right side is the same. In particular, the result is identical for the lines 1 and 2. Alternatively, we can use a residue 




(16) 




(17) 
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theorem. For the integration contour, we choose a path, which consists of the hne going from x — iy to x + iy (for 
instance, hne 2 in Fig. 1) and then clockwises along a semicircle, as again shown in Fig. 1. We then consider a limit 
y — >■ CXI, which ensures that S is infinitely small along the auxiliary semicircle, as follows from Eq. (1151) . The role 
of this semicircle is purely technical: it is needed to apply a residue theorem. The result of the integration depends 
only on how many poles of S{R) are enclosed by the contour. By enclosing the whole set of poles, we arrive to Z 
corresponding to the ground state of the initial quantum-mechanical problem, for which an equilibrium R is located 
to the left of the whole set of the one-electronic energies on a complex plane. Enclosing less number of poles leads to 
Z for excited states, for which an equilibrium R is known to be positioned between two one-electronic levels on real 
axis (quasi-contimium spectrum). Excited states are beyond the scope of the present paper and will be addressed 
elsewhere. Note that it is not necessary at all to consider paths consisting of straight lines (see, e.g., curve 3 in Fig. 

After integration, we readily obtain a very simple expression for the ground-state Z 

Z = cxp{-2eFo/V) £(-l)V" = exp{-2spjV) (1 - af^ . (18) 

n=0 ^ ^ 

By finding a logarithmic derivative of Z with respect to l/V, we recover a well-known expression of the single-pair 
energy 

Ei^2eFo-2n-^, (19) 
1 — a 

which also coincides with Eq. (1) for N — 1. 

This method yields the single-pair energy with error of the order of 1/p. The same however applies to the traditional 
method to derive Eq. (fT9|) . which relies on the replacement of the sum by the integral (see Section IV). 

The suggested scheme can be applied to the case of many pairs. Let us restrict ourselves for the moment to values 
of N such that N < Nq/2 for which the degree of the polynomial Yii j j<ii-^i ~ ^j)^ each Rj is smaller than that 
of the polynomial O^i UkC^^^ - Rj)- I* then follows from Eq. ^ that 5({i?j}), S{{Rj}) 'Z'j'^i Rj ^ at Im 
Rn — > ±oo for any i?„ from the set {Rj}. 

We consider a system of free charges in the equilibrium and then start to move this system as a whole out of 
equilibrium, in such a way that mutual distances between free charges stay the same, while the position of the center 
of mass R of this system changes. Now, we concentrate on the change 5W of W{{Rj}) upon this motion. By 
expanding W{{R,j}) in powers of the deviation 6R = R — R'^^^ of R from the equilibrium, we easily arrive to the 
equation 

oo 

SW ^Vj2'^n{SR)", (20) 



where 



(21) 



where i?^'''' are positions of free charges in equilibrium. Next, we assume that in the ground state all i?^"'' are located 
far enough from the line of fixed charges, so that distances between each free charge and this line are much larger 
than l/p. This assumption is rather natural for the thermodynamical limit. Besides, the method for the solution of 
Richardson equations developed earlier in Ref. f^] is also based on the same assumption, since it utilizes a continuous 
approximation for the positions of free charges. However, in contrast to f^], we do not introduce any additional 
requirements on the shape of a distribution of i?^"-* on the complex plane. 

In the limit — ?> oo, it follows directly from Eq. (|21|) that k„ depends on system volume as 

Kn ~ pNn^-'\ (22) 
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Moreover, K2 is purely real due to the mirror symmetry of the system of charges with respect to the real axis. Eq. 
(|22p implies that the absolute value of S is extremely strongly peaked near the equilibrium. Namely, if we go in the 
direction of the steepest descent of —Re W, then it is peaked within the interval of the width ~ If we now 

consider a ratio of integrals similar to the one in Eq. (jl2p . where the integration is performed for the position of the 
center of mass as well as for internal degrees of freedom of the system of charges, then we can find the equilibrium 
location of the center of mass, while the error will be in powers in as Eq. ([22]) suggests, so that it is negligible in 
the thermodynamical limit. We may formally extend the integration path for R from —ioo to +ioo, where S is zero. 
At the same time, the integrals entering the ratio are equivalent to the multidimensional integrals with respect to all 
i?'s. By using the Cauchy theorem, we may deform the integration path for each R, limits of integration being —ioo 
and +ioo. Thus, we are led to consider the ratio of integrals defined as in Eq. (jl2p . where an integration is performed 
for each Rj. The ratio E will allow us to reconstruct an information about the sum of all R's in an equilibrium, 
while this is the quantity of interest, since it is equal to the energy of the initial quantum- mechanical system. It is 
again convenient to introduce Z given by Eq. (jl4l) . where an integration is now performed for all variables Rj from 
—ioo to -|-ioo. The sum of equilibrium i?'s is also given by Eq. (fT3| . Eqs. (|13p and ([T4l) lead us to associate Z 
with the "partition function" for our classical system of charges, although such an analogy is certainly not complete. 
Nevertheless, we think that this term is far from being meaningless, in the view of Eq. A nontrivial feature is 

that the energy of the initial quantum mechanical system is determined by the logarithmical derivative of this classical 
" partition function" . 

Note that our approach has some analogies with thermodynamics. Indeed, for the system of many particles, it is 
often hopeless to resolve equations for equilibrium positions of each particle. However, such a detailed information is 
not necessary to understand many global properties of the system, so that a thermodynamical description becomes 
highly efficient. The major difference with the usual thermodynamics is in the fact that we are dealing with the 
unstable equilibrium and therefore the " partition function" is obtained by integrating only over half of the degrees of 
freedom and by using a complex- valued " probability" . We also would like to note that the method we here suggest 
turns out to have similarities with the large- iV expansion for the two-dimensional Dyson gas proposed recently in 
Ref. p^]. Perhaps, our approach may be also related to inverse problems in mathematics. For example, there exists 
the Radon transform, which allows one to reconstruct an unknown function by using integrals of this function along 
various paths This method is widely used in the tomography. 

Onc-dimcnsional integrals having denominator of the integrand of the form Jl£o('' ~ '^l^ilc integration for r 
is performed from —ioo to +ioo are called Norlund-Rice integrals. Here we deal with the multidimensional coupled 
integrals of Norlund-Rice type. At the same time, factor Yii j j<ii^i ~ ^j)^ the integrand makes them similar to 
Selberg integrals. Note that it is actually known that Norlund-Rice integrals can be transformed to binomial sums 
and vise versai^. 

In the next Section, we apply a technique based on binomial sums for finding the ground state energy for arbitrary 
N. However, before doing this, let us address one more point. 

We have imposed a restriction N < Nq/2. Configurations with N > Nq/2 can be handled by using a concept of 
holes, i.e., empty states in the potential layer. In the Appendix A, we show that the initial Hamiltonian is characterized 
by a duality between electrons and holes, which allows us to map states with N < 7Vo/2 to the states with A'' > Nq/2. 

III. GROUND STATE ENERGY THROUGH THE BINOMIAL SUM 

In this Section, we consider the ground state energy for N pairs, N < Nq/2. It is seen from Eq. (fTT|) that the 
"probability" S can be rewritten using a partial- fraction decomposition, in a manner, similar to the one used for the 
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one-pair problem, as 

S{{Rj = 2eFo - rj}) = exp 



V 



U3,l>3 



Nn 

E 

ni,n2...,njv— 



(23) 



It can be of interest to note that the sum given by Eq. (j23p may be rewritten in terms of the forward difference 
operators, as usual for binomial sums^i^. This means that Richardson equations can be represented in a similar way 
too. 

After substitution of Eq. (|23p to Eq. ([Ti)) and integrating in such a way that the integration path for each R avoids 
the whole set of poles, as for the one-pair problem, we obtain the binomial sum, given by 



Z = exp(-27VeFo/T^)z, 



where 



ni,Ti2...,njv— 



TV 

[](-i)"^-a" 



rii 



(24) 



(25) 



This sum again corresponds to the ground state of the quantum-mechanical system, for which none of equilibrium 
i?'s is located between two one-electronic levels in the real axis. 

Without the last factor in the RHS of Eq. (|25p , this sum reduces to the product of trivial binomial sums for each 
Uj, as the one of Eq. (jl8p . In order to tackle the sum with the coupling factor, we first provide some identities, which 
will be very useful in the derivation presented below. 

It is convenient to introduce Pochhammer symbol (or falling factorial) given by 



[n)q = n{n - l)...{n - g + 1), 



while (n)o = 1. It is then easy to obtain the following identity 

Nn 



/N \ 

z,,, ^ $](-l)"a" {nUNn - n), = (-l)V^(l - a)^--'^- 



{Nn-q-p)V 



(26) 



(27) 



where q + p < N^i. If we now consider a product of sums for each rij, every sum being similar to the one of Eq. (j27l) . 
we get 

Nn N 



-qi ,pi---^qN ,pn 



ni,n2...,njv— J — 1 



fNn 



N 



(28) 



iNn-qj-Pj)V 

Note that (i) only first two factors in the last line of Eq. ([28]) depend on the interaction constant V (through a); (ii) 
the dependence of both these factors on the two sets of numbers {qj} and {pj} is only through 9i a^nd ^J=i Pj^ 

which are just degrees of polynomials Y[j=i(j^j)qj ^^^^ 11^=1 (-^^^2 — "j)pjj respectively. These two observations turn 
out to be of a crucial importance. 

Our approach is to transform the initial coupled sum, as given by Eq. (|25p . into a linear combination of uncoupled 
sums, similar to the one given by Eq. (|28p . To make such an uncoupling, we note that, as known, Yii j lyji'^^i ~ '^i) 
can be rewritten as the determinant of the Vandermonde matrix 



n ("'-"j) = ^(K})=det 



1 


1 


1 


1 


m 


"-2 




■ nN 




9 


9 




n{ 




«3 ■ 


■ ri% 


N-l 








1 









(29) 
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Next step is to note that Yii j i>ji''^i ~ '^j) t)e also written as (— 1)^(^ ^''^^V{{Nq — rij}), since 
[{Nil ~ ni) — {Nfi — nj)]. Hence, we obtain the identity 

l,j,l>3 



ni - n, 



(30) 



Now we make use of the well-known rule that the determinant of the matrix does not change if we add a multiple 
of one row to another row. It is then easy to see that V{{nj}) can be rewritten in a "falling factorial" form as 



F(K}) = dct 



(»^i)o ('^2)0 ("-3)0 

(n-i)i (^12)1 ("s)! 

("-1)2 (^12)2 ("•3)2 

("Oat-I ('^2)Ar-l (?^3)^'- 



{nN)o 
{nN)i 
{nN)2 



N-l 



(31) 



while V{{Nci — nj}) can be represented in a similar form with rij changed into — nj. 

It is obvious that, using Eq. pip . V{{nj}) can be written as a linear combination of polynomials each having a 



form Y[j=ii''^j)qj ■ The crucial point is that for each term in this sum, 1j same: it is just equal to the 

degree of the polynomial V{{nj}). This degree is equal to the sum of degrees of polynomials from each row that is to 
-I- (TV— 1) = 7V(7V— l)/2. The same applies to V{{Na — nj}), which is represented as a linear combination of 



+ 1- 



polynomials of the form Y\j^i{Nn — nj)p. with Pj — N{N — l)/2. Now we see that Eq. (|3T|) together with the 

similar equation for V{{Nq — rij}) allows us to rewrite Yli j lyji^i ~ "-j)^ ^ linear combination of polynomials of 

the form (^Jlj^iC^Ogj) (^IljLil-^o — "-j)j3i^ with the same X^jLi Qj X^jLi-Pj each polynomial. At this stage, 
we immediately apply Eq. (|28p and obtain 



a^^(l-a)^^"-^(^-i)A(7V,7Va), 



(32) 



where A(N^ Nq) is some function of N and A^o, which is irrelevant for the determination of the quantum- mechanical 
energy, since the later is given by the logarithmic derivative of Z with respect to 1/V. Finally, by finding this 
logarithmic derivative and by taking into account Eq. (|24p . we easily arrive to Eq. (1) for the ground state energy. 

Eq. ([5^ can be derived using more formal way of writing. We first note that V{{nj}) can be expressed in the 
following form 



Vi{n,}) = det{V,,p) = det[(nj,),_i] 



AT 



ejii2 . . .jN ('^ii ) ('^i2 ) 1 ■ • ■ ) Af- 1 > 



(33) 



where s 



jij2-..jN is the Levi-Civita symbol. Similarly, we can write 

N 



N 



jl ,31 , — ,jN = l il J'2 ' ■ ■ ■ J'jv = 1 



(34) 



By definition, the Levi-Civita symbol is nonzero only for the set of its indices all different. This means that, for 
nonzero terms in the double sum in the RHS of Eq. ([34]), there should be no repetitions in the set ii, J2, ■■■tJn and 
also in another set jj, j^, j^v- ^^^^ ^'i- and see that each nonzero term of the double sum gives the 
same dependence on V, when substituted into Eq. ([M)) . Thus, we again arrive to Eq. (15^ . 

In order to reach some qualitative understanding of the result, we have obtained, let us come back to the expression 
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of " probability" , as given by Eq. (fTTj) . We rewrite it as 



S{{R, = 2ep, - r,}) = exp(-27Ve;^„/y) 



1 



exp 



V 



det 



1 1 ... 1 ■ 

nLo('^i + f) it=,{r2 + '-f) ... n°=o('^A^ + f) 



det 



1 1 



2n^ 
P ^ 



(35) 



It is seen from the above equation that the "probabihty" S{{Rj}) can be represented in a factorized form as a sum 
of products of N "probabihties" , each of them being the one for a single pair. Each single pair is however placed into 
its own environment with a band of one-electronic states removed from the bottom of the potential layer and another 
band of states removed from the top of the potential layer, as a result of the Pauli exclusion principle. The energy 
of a single pair in the state with n levels removed from the bottom of the potential layer and m states removed from 
the top of the layer is obtained from Eq. by substitution ep^ + n/ p, il ^ — {n + m)/p. Hence, it follows 

from Eq. ((35)) that the sum of single-pair energies is the same for each term of the factorized " probability" , although 
sets of single-pair energies for different terms are not identical. 

We therefore can think that the original system of N pairs, feeling each other through the Pauli exclusion principle, 
is a superposition of "states" of N single pairs, but each pair is placed into its own environment with bands of the 
states deleted both from the bottom and from the top of the potential layer. Moreover, the sum of energies of these 
single pairs for each "state" of the superposition is the same, which is, probably, a consequence of the constant density 
of states. This understanding provides an additional nontrivial link between the one-pair problem, solved by Cooper, 
and many-pair BCS theory. 

Note that we see here a quite close analogy with the well-known Hubbard- Stratonovich transformation which 
enables one to represent the partition function for the system of interacting particles through the partition function 
for the system of noninteracting (single) particles, but in the fluctuating field. We also would like to mention that 
some of the "probabilities" appearing in Eq. ([55]) can be negative; one should not be confused by this fact, since 
negative "probabilities" are known to appear in problems involving fermions. In particular, they lead to the well- 
known negative sign problem arising when trying to apply Monte Carlo methods to compute partition functions. 
This problem becomes severe in the thermodynamical limit. Actually, we here see some reminiscence of the same 
problem. Namely, if we try to integrate factorized S through paths crossing individual saddle points of each single 
pair, we immediately see that the absolute value of resulting Z is going to be much smaller than the result of the 
integration for each term due to the " determinantal" structure of Eq. (|35l) . The easiest way to see it is to consider the 
two-pair problem. The situation is getting more and more difficult with the increase of N. Fortunately, by applying 
the method based on binomial sums, we can circumvent all these difficulties and to evaluate Z exactly. 



IV. GROUND STATE ENERGY THROUGH THE NORLUND-RICE INTEGRAL 



In this Section, we present another method, which allows us to find the ground state energy as an expansion in 
powers of pairs density. In practise, this method is much less powerful compared to the one described in the previous 
Section. We here calculate just few first terms of the expansion, since calculations become more and more heavy with 
increasing number of terms. One of the main aims of this Section is actually to establish a link between our approach 
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and the dilute-limit approach of Refs. fiiii^]. The method presented in this Section turns out to be closely related to 
the method of Refs. fii^]. 

The general idea is not to transform the initial multidimensional integral of Norlund-Rice type into a binomial sum, 
but to tackle it by using a saddle-point method. For z, we have 



z = I dr — -J — — — exp ' 



2m V 



V 



(36) 



where 



dr — I dri ... / dr 



FN- 



Quantum- mechanical energy E is expressed through the logarithmic derivative of z with respect to 1/V as 

N 



E = 2NeFa 



dln{z) 

WJv) 



2Nep^ --idr\Y: 



^ / n n " { 



^1=1 



(37) 



The factor 0; j j<i(^i ~ ^j)^ is responsible for the coupling between r^-'s. In order to see the effect of this coupling, 
we use the path of integration for each rj , which passes through the stationary point corresponding to a single pair 
(line 1 in Fig. 1). This point yields the one-pair energy, as found by Cooper. 

Because of the symmetry between all r^'s in the RHS of Eq. ((37)) . we may replace X^jLi ''j by Nri. Then, E can 
be rewritten as 



E = 2NeF,-N^, 



(38) 



where 



drr 



„ Y\i.j.j<i{n~rjf 



2^1=1 



V 



drr'^ 



N / Nn 



2mj 
P 



rrij—O 

The position of the saddle point is given by the solution of the following equation: 



(39) 



-fir) ^ - 
dr dr 



r 
V 



El- 



2m 
P 



V ^ r + — 

m=0 ^ P 



0. 



(40) 



We replace sum by the integral that is justified in the large-sample limit, and find the solution of Eq. (|40| . r — tc, 
which gives the binding energy of an isolated Cooper pair 



2rj 



1-(t' 



(41) 



which actually has already been found in Eq. (|19p through the binomial sum. We now introduce a new variable t 
defined as it — r — Ec- Next, we expand /(r) in a Taylor series around r = ec and rewrite it as 



/(e, + ^^) = ./(ec) +¥'(*), 



(42) 



where tf(t) reads as 



(43) 
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Derivatives (cc) can be easily found from Eqs. (|40l) and (|411) as 



/(")(Oc.(-ir(m-2)!^ 



p / 1 — cr 



2 V 217ct 



(1 - (7™-^). 



(44) 



Notice that we again have replaced sums by integrals, when deriving Eq. (j44p . 

It is seen from Eq. (|44l) that /'•^-'(ec) is positive, since a < 1. Thus, the direction of the steepest descent for 
exp(/(r)) is perpendicular to the real axis in r complex plane, as expected (line 1 in Fig. 1). Hence, we assume that 
t changes from — oo to +oo. Next, we introduce a rescaled variable x instead of t 



Within new notations, / reads as 



where 



(45) 



(46) 



m=2 



m(m — 1) 1 — cr 

where is a small dimensionless constant which is inversely proportional to the square root of the system volume 

1 



(47) 



(48) 



Note that a2 = -1/2. 

The existence of the small constant S is crucial in our derivation, since we are going to calculate the ground state 
energy as an asymptotic expansion in powers of d. We will however see that 5 is going to be coupled with the number 
of pairs, so the resulting expansion of the energy density in the thermodynamical limit will be in powers of pair 
density, which is proportional to NS^ ~ N/p and is no longer small. Such an unusual feature is due to the fact that 
we here deal with the multidimensional integral, while integration is performed through the single-pair saddle point. 

Next, we rewrite Eq. (|38l) as 

2na $1 



E ^ 2NeFo - - iN6 



1 - cr $Q 



where 



dxxl 



N 

J|exp 



(49) 



(50) 



where integration for each Xj is performed from — oo to +00. 

Let us concentrate on $1. It can be integrated by parts with respect to Xi in the following manner: 



$1 = I dx 



n (x,-a;,fe^^-"'"-" 



(51) 



The derivative in the integrand in the RHS of the above equation reads 



+6 amO Xi 



m=2 



n (x,-a:,fe^^-"'"-" 



(52) 
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where 



( ,■^m+l 1 _ 

.i(m + i)r(™-i) = i^^^ 



^m+U'"- ±yu ■ -— . (53) 

m 1 — (7 

The first term in Eq. ([5^ vanishes identically after we substitute it into Eq. ([?T|) and perform integration over all 
Xj's, due to the symmetry reasons, which can be easily verified upon the mutual exchange xi i — > Xj-^. The second 
term, after the integration, yields an expansion of $i through $2, '&3, etc. with coefficients, which represent higher 
and higher powers of 6 

00 

$1=5^ a„,(5™-2^™. (54) 

We now focus on the first term in the expansion from the RHS of Eq. (|54p . which is linear in S. $2 can be handled 
in the same manner, as the one used for $1. Namely, we make integration by parts with respect to xi. We then 
obtain an equation, similar to Eq. (j5ip . except that the expression in square brackets must be multiplied by xi. The 
derivative of this expression can be written as 

n (x,-x,)'e^"-"'"^" 
+2x,J2{x,-x,,) n {xi-x.fe^^-^'^-^^ 

0\07^(ljl) 

+<5 ^ a™5™-2<+i J] (xz-^,)2e^--3"™-r. (55) 

m=2 l,3,]<l 

After integration, the first term in the RHS of the above equation gives $o- The second term is the sum of — 1 
terms. By performing a mutual exchange xi < — > Xj^ , it is easy to see that each of these terms also gives $0 after the 
integration. The third term leads to the linear combination of $m's with m starting from 3. As a result, $2 reads 

00 

<^>2=N<^>o + 6j2 a™<5"-2$,„+i. (56) 

m=2 

At this stage, we substitute Eq. ()56p to the expression of $1, as given by Eq. ([M)) . and obtain 

00 

m=2 

= SN<i>aa2 + 0{6'^), (57) 
where a^^j^ is a "new" expansion coefficient, given by 

(58) 

Eq. ()57p provides $i/$o in lowest order in S. It is then substituted to the expression for the energy, given by Eq. 
. Together with Eq. ([551) for m = 2 and the definition of S it leads to 



E = E, + ^\±^ + 0(p-^^A, (59) 

p 1 — (7 V / 



where En is given by Eq. (1). The second term in the RHS of Eq. ()59p is underextensive, so that it can be dropped. 
Of course, one has to keep in mind that the third term, O (p^'^/^), is not necessarily small due to the coupling with 
N. 

Thus, we have calculated the correction to the energy of A'^ noninteracting pairs in the lowest order in S, which 
gives nonzero contribution, i.e., in 6"^. This correction turns out to be proportional to N{NS'^) (due to the coupling 
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with N) that is extensive. Moreover, this first correction already gives the total energy within all relevant (extensive) 
terms, as evident from the results of Section III. The task is now to prove that the third term in the RHS of Eq. (|59|) . 
produces underextensive contribution only. This problem is addressed in the Appendix B, where we show, through 
rather tedious calculations, that two next terms in S are indeed underextensive. We stress that we are unable to 
present a complete proof involving all terms of the expansion. The underextensivity we have revealed follows from 
"magic rules" for coefficients Um, which are of course linked to the constant density of states. Thus, the method 
based on calculation of the binomial sum, which enabled us to derive the expression for the ground state energy 
directly within all extensive terms in a simple way (nonperturbatively in S), appears to be much more efficient than 
the method of the present Section. 

We also wish to stress that the method of saddle-point evaluation of Norlund-Ricc integral through the single-pair 
solution is quite similar to the method of analytical solution of Richardson equations presented in Ref. [^^] and 
advanced further in Ref. \^]. In Ref. fi^], all Richardson energy-like quantities i?'s have been expanded around a 
single-pair solution in seria involving small dimensionless parameter analogous to 6. By keeping only the lowest-order 
contribution, Eq. (1) has been derived. The lowest-order approximation a priori makes this method restricted to 
the dilute regime of pairs only, when their density is small. However, in Ref. p^] it was demonstrated that next 
term in the small dimensionless parameter vanishes identically due to the first " magic cancellation" , which implies 
that the dilute-limit result of Ref. fi^] is likely to be universal. We have arrived to the similar conclusion within our 
framework, but we have also revealed the second "magic cancellation" which exists for the next-order term (see the 
Appendix B). 

The method of integration through the single-pair saddle point, as presented in this Section, has some disadvantages 
(in addition to its obvious technical complexity compared to the method based on manipulations with binomial sums) . 
Namely, it relies on replacing sums by integrals, as we did to derive Eq. (|1T|) . More rigorously, sums of this kind 
should be expressed through F-functions. By keeping a leading order term in the expansion of these F-functions in 
1//9, one obtains the approximation, used here. Dropping all the remaining terms introduces errors of the order of 
5^ ~ 1/p in many places. Fortunately, this approximation does not lead to any pathologies for extensive quantities, 
but it produces artificial underextensive terms (not considered here). In addition, the effect of the reduction of a pair 
binding energy in a many-pair configurationi^ compared to the isolated pair is clearly of discrete origin, while, within 
this approach, it is recovered through the continuous approximation, i.e., in a very indirect way. 

V. CONCLUSIONS 

Within the exact mapping of Richardson equations to the system of interacting charges in the plane, we suggested 
to switch from the "energy" of this system to the "probability" for charges to occupy certain states in configurational 
space at the effective temperature given by the interaction constant. The effective temperature thus goes to zero 
when system volume goes to infinity. We introduced a " partition function" , from which the ground state energy of 
the initial quantum-mechanical many-body problem can be found. This approach leads to the emergency of quite 
rich mathematical structure. The "partition function" has a form of a multidimensional integral, similar to Selberg 
integrals. For the model with constant density of energy states, the structure of the integrand implies that it can 
be also considered as an integral of Norlund-Rice type. The most efficient way to evaluate it is by transforming the 
integral into a binomial sum, where the coupling between variables is due to the determinant of the Vandermonde 
matrix. Using properties of this matrix, we managed to evaluate the "partition function" exactly in a rather simple 
way and to find the ground state energy, which is described by a single expression all over from the dilute to the dense 
regime of pairs. This expression does coincides with the mean-field BCS result. 

We also provided a qualitative understanding of the obtained result in terms of "probabilities". Namely, the 
"probability" of the system of N pairs, feeling each other through the Pauli exclusion principle, to be in a certain 
state can be represented in a factorized form as a linear combination of terms. Each of them is given by the product 
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of "probabilities" for N single pairs, but every single pair is placed into its own environment, which is identical to 
the initial one, except that two bands of the one-electronic energy states - both from the bottom and from the top of 
energy interval - are absent due to the Pauli exclusion principle. Moreover, we find that although these environments 
are different, the sum of energies of the single pairs is the same for all terms of factorized "probability". 

Finally, we presented another method for evaluation of the Norlund-Rice integral by integration through the saddle 
point corresponding to a single pair. This method turns out to be much less efficient; it exactly corresponds to the 
dilute-limit approach for the solution of Richardson equation proposed very recently in RefsJ^ii^. Nevertheless, we 
managed to calculate several initial terms of the expansion of energy in pairs density. 

The suggested method is rather general and we believe that it can be applied to other integrable pairing Hamilto- 
nians, which have electrostatic analogiea^iiS^, as well as to situations with nonconstant density of energy states. 
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Appendix A: Electrons-holes duality 

We introduce creation operators for holes as b\^^ — Okf and h\^^ — with corresponding rules for destruction 
operators. These holes should not be confused with holes in the Fermi sea of normal electrons, as appear in BCS 
theory. By using commutation relations for fermionic operators, it is rather easy to rewrite the initial Hamiltonian in 
terms of holes as 

= ^ 1 + 2 ^ ek - E (^k - (&U^kt + hi,h^,) -VY, &Lt^-k'i&-k4&kf (Al) 

k k k k,k' 

The first two terms of the RHS of Eq. (jAl[) are scalars. Moreover, the first term, V 1, is inderextensive and can be 
dropped. The fourth term in the RHS of Eq. (|A1|) is fully identical to the interaction potential in terms of electrons, 
given by Eq. To better understand the role of the third term, we introduce defined as = + — Ck- 

In the case of a constant density of states, is 0, 1/ p, 2/p, VI: just counts states starting from the top of 
potential layer towards its bottom. Hence, — (ck — V) can be rewritten as — {e Fq -\- ^ — V) . A similar term in the 
Hamiltonian for electrons, given by Eq. ([3]), contains factor Ck +£Fo- Thus, we clearly see a duality between electrons 
and holes: The ground state of N pairs is equivalent to the ground state of A^o — N holes, while the energy of the 
latter is equal to the energy of Ao — A^ pairs, with changed into — (£_Fo + fi — V), plus the kinetic energy of the 
potential layer completely filled, given by 2^j.ek. This means that we can avoid considering configurations with 
A' > A'o/2 by switching to holes. Cooper pairs in this case being constructed out of them. Moreover, one can directly 
check that the expression of the ground state energy given by Eq. (1), which we are going to prove, satisfies, within 
underextensive terms, the above duality criterion. Therefore, it is sufficient to consider the case A^ < Nq/2 only. The 
existence of the duality between electrons and holes in the energy spectrum has already been noted in Ref. [^^'^^], 
although its origin, at the level of the Hamiltonian, stayed unclear. 

Appendix B: Cancellation of higher-order terms in pairs density 

We focus on $3 and rewrite it through the integration by parts, as has been done for $1 and $2 

00 

$3 = 2Ar$i + 6 am'5™"'*™+2. 
Next, we substitute Eq. (|Bip to the first line of Eq. (|57l) and get the following expression of $i/$o 

*1 TAr 1 f-, , ^ (2) cm-3*m+l\ 

— = SNa2 TTv 1 + -r^ > al'^S'"- ^—r^ , 

*o 1 - 2Smai'^ V ^3 $0 J 



(Bl) 



(B2) 
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where 



(2) _ ^(1) , ^(1 



^m+1 ^ 



'_^^ + a);'a,n-i- (B3) 



We now consider the first term of the sum in the RHS of Eq. (|B2I) . Again integrating by parts, we obtain 

oo 

$4 ^ 27V$2 + iV$i,i + S J2 am-2<5"-^$™+i, (B4) 

m— 4 

while $m,Ti is given by Eq. ([SO)) with replaced by x'^X2- ^i,i is evaluated as 

oo 

= -^o + SY^ a,„+i,5™-i$i,^+i. (B5) 

m—l 

Now we substitute Eq. (jB5l) for and Eq. ((56|) for $2 to Eq. (|B4[) for $4, and the resulting equation for $4 - to 
Eq. ()B2[) for $1. The obtained expression of $1 can be written as a sum of two contributions 



$i=Gi+G2, (B6) 

where 

G, = 6N^^^ ( 1 + ^S'n] , (B7) 

1 - 2a^^^>d^N \ a2 J 



r4 / 00 00 \ 

1 - 2ai'^Sm V^4 ^2 J 

where 

a^+i - a„_,_i + 04 am-2- (B9j 

It is can be checked that 

af^+ 4^^02 = (BIO) 

for any a. Therefore, due to the exact cancellation between the numerator and denominator, Eq. (|B7[) for Gi is 
reduced to 

Gi = {SN)^oa2- (Bll) 

This means that terms proportional to N6{N6'^) — which are present in Gi and not in G2, are absent in the 

expansion of $i/$o in powers of 6. As clearly seen from Eq. (P^. these are exactly the terms, which give contribution 
of the order of N{N/ p)"^ ~ N{N6'^)'^ to the energy, or equivalently, terms in square of pairs density {N / pY) to the 
energy density. Eq. (|B10I) constitutes first "magic cancellation" rule and it is fully equivalent to Eq. (68) of Ref. fi^]. 

Let us now make some comments. In general, at each step of our procedure, we express through the linear 
combination of all possible terms of the form '^mx,m2,...^ with the sum of nonnegative integer m's equal to m — 2 and 
prefactors, independent on 5, plus expansion which involves (5"$m+„, n starting from 1. We then have to proceed 
with all these "new" <I'mi,T?i2....'s to lower the sum of m's by 2 at each iteration and finally to bring them either to 
$0 or to $1, depending on a parity of m. And so on. Of course, this recursive tree-like procedure becomes more and 
more tedious when m is increasing. That is why we were able to trace only few first terms of the expansion. 

Now we provide a sketch for the derivation of the next-order correction to the energy, which will be of the order 
of N{N/ pY. For that, we first have to consider $5. We express $5 through $3 (given by Eq. (jBl[) ). $1.2, and $„ 
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for n's starting from 6. $1^2 is expressed through $1 and $2,n, with n starting from 2. Then, we perform similar 
manipulations for $6- Next, we substitute all these quantities to Eqs. (|B6P - (|B8|) for $1 and obtain more sophisticated 
fraction than the one in the RHS of Eq. (|B2[) which now includes higher powers of 5. 

By expanding this "new" fraction in 5, we finally ensure that terms of $i/$o in N5{N5'^Y — vanish due to 

the second "magic cancellation" rule, given by 



(B12) 



where 



a. 




(B13) 



Again, Eq. (|B12|) is fulfilled for any <t. Then, it follows from Eq. (j49l) that terms in N{N/p)^ ^ N{N6'^)^ vanish in 
the expression of the quantum-mechanical energy. 



